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INTRODUCTION 


This  document  constitutes  the  final  report  of  the  second  phase  of 
an  investigation  into  convection  in  the  marine  atmosphere  and  into  topics 
pertaining  to  radar  sensing  of  the  atmosphere.  In  a  short  introductory 
phase  reported  in  the  contract  final  report  (Freeman  1979),  we  reviewed 
relevant  theoretical  models  of  marine  convection,  surveyed  the  climato¬ 
logical  frequency  of  convection  in  the  North  Atlantic,  and  offered  recom¬ 
mendations  on  the  modelling  of  triggered  convection. 

The  current  investigation  has  also  had  a  limited  scope  of  study; 
the  principal  objective  is  to  develop  a  technique  for  evaluating  the 
conditions  under  which  boundary  layers  are  unstable  and  for  computing  the 
rate  of  growth  of  unstable  motions.  In  recognition  of  the  limited  time 
and  funds  available  for  this  study  we  restricted  our  considerations  to 
simple  models  of  the  boundary  layer,  which  can  be  evaluated  without  large 
computer  expense  and  complication.  Nevertheless,  it  has  been  possible  to 
include  in  this  framework  all  of  the  physical  effects  which  can  contribute 
to  the  perturbed  motion  of  the  marine  boundary  layer.  We  take  into  account 
the  influence  on  linear  stability  of  wind  shear,  static  stability,  Coriolis 
force,  turbulent  diffusion,  large-scale  convergence,  and  the  geostrophic 
wind.  A  procedure  has  been  developed  for  solving  the  linear  perturbation 
equations  to  obtain  the  complex  eigenvalues  and  eigenfunctions  correspond¬ 
ing  to  a  given  unperturbed  state  of  the  boundary  layer.  The  equations  are 
sufficiently  general  that  realistic  profiles  of  wind  speed,  density,  and 
diffusivity  can  be  accomodated.  We  have  provided  alternative  sources  of 
the  profiles  of  mean  quantities;  they  can  be  prescribed  as  analytic 
functions,  or  they  can  be  supplied  as  numerical  tables  derived  either  from 
data  or  from  a  boundary  layer  model. 


In  our  report  to  follow  we  describe  this  new  computational 
technique  and  the  results  obtained  with  it.  We  also  report  on  related 
calculations  with  the  boundary  layer  computer  code  SIGMET,  some  studies 
of  radar  properties  of  the  marine  boundary  layer,  including  ducting,  and 
an  investigation  of  atmospheric  internal  wave  propagation. 

The  report  is  divided  into  four  sections.  In  Section  A  we  consider 
the  mean  structure  of  the  boundary  layer  as  displayed  by  the  SIGMET 
computer  code;  stability  of  the  boundary  layer  is  discussed  in  Section  B, 
where  formulation,  mathematical  implementation,  and  numerical  studies  are 
presented.  This  section  contains  the  major  development  of  the  study.  In 
Section  C,  radar  properties  of  the  marine  atmosphere  are  reviewed  and  in 
Section  D  a  study  of  atmospheric  internal  wave  propagation  is  reported. 

We  believe  it  should  be  possible  to  associate  radar  surface  ducting  and 
internal  wave  activity  with  other  properties  of  the  marine  atmosphere. 


A.  MEAN  STRUCTURE  OF  THE  MARINE  BOUNDARY  LAYER 

One  of  the  capabilities  which  we  require  for  the  analysis  of  the  boundary 
layer  is  a  description  of  the  profiles  of  temperature,  humidity,  wind  field, 
and  diffusivity.  These  quantities  evolve  through  diurnal  and  inertial  cycles 
in  response  to  external  forcing  parameters,  such  as  sea  surface  temperature, 
geostrophic  wind,  large-scale  convergence  or  divergence  of  the  horizontal 
wind,  and  short-  and  long-wave  radiation.  They  are  also  influenced  by  initial 
conditions  of  atmospheric  temperature,  humidity,  and  wind  field.  In  response 
to  these  quantities  the  boundary  layer  changes  with  time  as  determined 
by  the  Navier-Stokes  equations  and  an  approximation  (based  on  second-order 
closure)  to  the  fields  of  several  second-order  variance  and  covariance 
quantities  formed  from  turbulent  fluctuations. 

The  above-described  model  determines  the  laminar  boundary  layer  in 
terms  of  the  average  or  unperturbed  profiles  of  the  atmosphere.  Deviations 
from  the  mean  motion  have  been  parameterized  in  terms  of  the  turbulence 
description  as  they  affect  the  mean  motion  itself.  This  model  is  required 
as  a  starting  point  for  a  more  detailed  model  (  see  Section  B)  which  describes 
the  response  of  the  laminar  layer  to  small  perturbations.  It  is  also  the 
starting  point  for  the  calculation  of  derived  profile  quantities  (see  Section 
C)  which  characterize  the  propagation  and  scattering  of  microwaves  in 
the  atmospheric  boundary  layer. 

The  unperturbed  boundary  layer  is  described  by  the  appropriately 
averaged  Navier-Stokes  equations.  These  are  given  in  Eq.  B-2  and  in  more 
detail  in  an  appendix  to  our  proposal  to  the  Naval  Research  Laboratory 
(Freeman,  1978).  In  the  latter  report  the  formulation  of  the  SIGMET 
computer  code  embodying  these  equations  is  also  described  and  several  examples 
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of  calculations  with  the  code  are  presented. 


In  support  of  the  current  investigation  we  have  carried  out  several 
additional  calculations  of  the  marine  boundary  layer  with  SIGMET.  These 
calculations  have  been  used  as  input  data  for  the  stability  calculations  of 
Section  B  and  to  provide  the  mean  profiles  from  which  the  microwave  modified 
refraction  coefficient  and  structure  function  were  calculated  in  Section  C. 
The  initial  data  for  the  calculation  are  shown  in  Fig.  A-l,  and  the  profiles 
of  velocity,  humidity  and  potential  temperature  at  a  time  (0108  local  time) 
eighteen  hours  after  the  initial  time  are  presented  in  Fla.  A-2  through  A-5. 

Finally,  we  have  incorporated  the  terms  describing  the  large-scale 
convergence  or  divergence  of  the  horizontal  winds  into  SIGMET.  This  term 
can  be  prescribed  as  an  arbitrary  function  of  altitude.  From  it  the  vertical 
velocity  of  the  air  column  is  calculated,  which,  in  turn,  determines  the 
vertical  advection  terms.  Applications  of  this  version  of  SIGMET  will  be 
made  at  a  later  time. 
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POTENTIAL  TEMPERATURE 
SPECIFIC  HUMIDITY 


FIGURE  A-l.  Initial  potential  temperature  and  specify  humidity  profiles 
used  ^or  a  SIGME"  atmospheric  boundary  layer  calculation. 
Local  time  s  OIOS  hours. 


B.  STABILITY  OF  BOUNDARY  LAYER  FLOW 

In  reviewing  literature  citations  to  determine  what  physical  and 
mathematical  approximations  have  been  investigated,  we  found  several 
studies  of  free  convection  {without  wind  shear)  and  others  of  shear 
instability  in  a  neutrally  stable  atmosphere.  The  dependences  of 
diffusivity  and  wind  on  altitude  were  highly  idealized,  and  the  role  of 
a  time-dependent  unperturbed  state  had  received  little  attention. 

We  have  formulated  the  perturbation  problem  more  generally  to  take 
account  of  these  simpl ications.  The  perturbation  equations  are  consistent 
with  those  solved  by  the  computer  code  SIGMET  which  evaluates  the  evolution 
of  profiles  of  wind,  temperature,  humidity,  etc.,  in  a  turbulent  boundary 
layer.  These  perturbations  may  be  unstable;  the  degree  of  instability 
depends  in  general  on  external  parameters,  such  as  heat  flux  and  geo- 
strophic  wind,  and  on  the  wavelength  and  direction  of  the  perturbation. 

For  each  such  combination  a  number  of  vertical  mpdes  are  permissible. 

The  initial  mathematical  and  numerical  formulation  was  chosen  to 
facilitate  checking  with  previously  obtained  results  and  to  permit  assess¬ 
ment  of  the  accuracy  of  the  numerical  solution.  We  assume  for  this 
purpose  that  the  unperturbed  state  is  steady  and  that  the  parameters 
derived  from  the  unperturbed  solution  are  provided  on  a  finite  grid  of 
representative  altitudes,  as  is  the  case  using  SIGMET.  These  points  are 
not  necessarily  equally  spaced.  Since  the  coefficients  are  independent  of 
time,  we  may  assume  that  the  time  dependence  of  the  perturbation  amplitude 
corresponds  to  a  constant  phase  speed.  If  the  phase  speed  is  real,  the 
perturbation  is  a  freely  traveling  wave  (for  example,  an  internal  wave), 
but  if  there  is  an  imaginary  component,  the  wave  is  either  damped  or 
amplified.  The  evaluation  of  the  complex  phase  speed  (or  frequency) 


10 


constitutes  the  basic  problem  of  determining  the  linear  perturbation 
solution.  In  the  above  formulation  the  solution  is  expressible  as  a 
matrix  eigenvalue  problem  containing  complex  coefficients.  The  number 
of  modes  is  determined  by  the  number  of  altitude  intervals;  the  higher 
modes  of  the  continuous  problem  are  suppressed  by  the  finite  difference 
representation.  The  eigenfunctions  corresponding  to  a  given  eigenvalue 
provide  the  altitude  dependence  of  the  linear  perturbation. 

In  the  following  sections  the  formulation  of  the  incompressible, 
steady  state  linear  stability  problem  is  described.  The  calculation  takes 
into  account  the  effects  of  wind  shear,  atmospheric  stability,  large-scale 
pressure  gradient  and  convergence,  the  Coriolis  force  and  turbulent 
diffusion;  these  are  included  in  the  equations  given  in  Section  B-l.  The 
differential  equations  contain  derivatives  of  order  four  and  lower;  in 
Section  B-2  these  derivatives  are  approximated  by  difference  equations  of 
order  four.  The  SIGMET  spatial  grid  is  not  uniform  in  vertical  spacing 
in  order  to  provide  higher  resolution  near  the  surface  where  the  wind 
shear  is  large;  we  account  for  this  nonuniformity  by  employing  a  coordinate 
transformation  described  in  Section  B-3.  Boundary  conditions  and  methods 
of  solution  of  the  complex  matrix  eigenvalue  equations  are  presented  in 
Sections  6-4  and  B-5.  The  resulting  computer  code  (PERT)  has  been  used  to 
evaluate  the  stability  of  a  large  number  of  representative  boundary  layer 
flows.  To  test  the  correctness  and  accuracy  of  the  solutions,  comparisons 
were  made  with  several  reported  stability  analyses;  these  are  discussed 
in  Section  B-6.  A  large  number  of  parameters  governing  stability  of 
atmospheric  flows  can  be  explored  with  the  PERT  computer  code.  In 
Section  B-7  some  of  these  have  been  investigated,  including  the  simultaneous 
effects  of  static  stability,  wind  shear,  Coriolis  force,  and  large-scale 
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B. 1  The  Equations 

The  Navier-Stokes  equations  for  a r  incompressible  stratified  fluid 
take  into  account  the  conservation  of  momentum  in  a  viscous  fluid  and  the 
equation  for  the  conservation  of  mass  in  the  presence  of  diffusion.  These 
equations  also  include  the  effect  of  the  Coriolis  force  in  addition  to 
the  pressure  in  the  fluid.  In  a  fixed  coordinate  system  x  (positive 
eastward),  y  (positive  northward),  and  z  (positive  vertically  upward)  the 
equations  are 


H  g.  uju  +  v|2  +  w|£  -  fV  +  f'W  ♦  1  I7  *  T  •  (K-U), 

ct  cX  dy  cZ  k  cX 


3V  .  ,.sV  9V  .  „3V  .  .  1  IF  .  -  , 

It  +  L,55  *  v37  +  "sT  +  ™  +  c  3y  '  ' 

3W  4.  ti3W  4.  v3W  4.  .  IMI  .  I  H  +  o  S  "  • 

3t  +  u^  +  v37  +  wsI  fU  ~T7  +  g  • 


H+Uft  +  Vf7  +  Wtfs7  •  (K'7M' 


SU  +  1Y  +  3W 
lx  5y  +  3z 


B- ( 1 ) 


Equation  B- ( 1 )  contains  the  wind  velocity  components  (U,  V,  W),  the 

pressure  P,  the  density  p,  the  Coriolis  parameters  f  =  23sinc,  and 

f‘  =  2Qcosp,  and  the  (primarily  turbulent)  viscosity  K  and  diffusivity  K'. 

Our  approach  is  to  consider  solutions  to  B- ( 1 )  which  represent 
horizontally  stratified  boundary  layer  flow  (the  unperturbed  solution) 
and  to  investigate  its  linear  stability  against  perturbations  involving 
horizontal  inhomogeneity.  The  equations  for  the  perturbations  are 
linearized  by  neglecting  all  nonlinear  terms. 

The  mean  flow  equations  neglect  all  horizontal  variations  except 
for  an  assumed  large-scale  (synoptic)  pressure  gradient  and  convergence 
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or  divergence  of  the  horizontal  wind  speed.  These  ter-s,  which  can  be 
functions  of  altitude,  are  assumed  to  be  derivable  fro-’  synoptic  weather 
data.  Including  these  terms,  the  mean  flow  equations  are 


li +  wl¥  -  f  (v-v  +  f'“  *  I?  <k  f). 


|V  +  W|Y  +  f  ,u-Ug)  .  4  (K  | |) 


92 


||  +  -  f'U  +  i|£  + 


Tt  T  "7z 


3  SW. 
9z  (K  9z  ' 


+  W  =  A  <K*  t#>. 


T? 


9  z 


9z' 


9W 

7*  =  "  °* 


B-(2) 


In  the  third  equation  for  the  vertical  component  of  momentum  it  is 

permissible  to  neglect  the  terms  containing  the  velocity  because  over  a 

flat  surface  they  are  very  small  compared  with  the  acceleration  of 

gravity  g,  and  the  term  containing  the  hydrostatic  pressure  gradient 

.  In  this  (hydrostatic)  approximation  the  vertical  momentum  equation 
3z 

becomes 


The  final  equation  gives  the  vertical  wind  speed  by  imposing  the 

surface  boundary  condition  ^  *  0  and  integrating  the  horizontal  divergence 

0  *  |H  +  ,  which  is  a  known  function  of  altitude.  Other  quantities 

3x  dy 

appearing  in  B-(2)  are  the  geostrophic  wind  components 
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assumed  to  be  known  functions  of  altitude.  The  turbulent  viscosity  K  and 
diffusivity  K‘  are  determined  by  appeal  to  turbulence  theory.  In  1-D  SIGMET, 
which  solves  the  above  system  of  equations,  a  set  of  auxiliary  equations 
is  solved  to  determine  the  turbulent  diffusivities.  These  are  obtained  by 
second  order  closure  hypotheses;  a  prognostic  equation  for  the  turbulent 
kinetic  energy  and  algebraic  equations  for  derived  turbulent  quantities 
result.  In  the  following  applications  we  are  able  to  use  the  1-D  SIGMET 
code  to  provide  data  on  the  mean  flow  solution,  including  the  turbulent 
diffusivities.  In  doing  so,  we  require  a  number  of  quantities  (such  as 
f,  D,  Ug,  Vg,  surface  boundary  conditions,  etc.)  which  specify  the  problem 
in  terms  of  local  and  synotic  parameters.  Their  values  determine  whether 
the  unperturbed  atmosphere  is  stable  or  unstable  and  the  amount  of  wind 
shear  and  turning  in  the  boundary  layer.  A  more  detailed  description  of 
the  SIGMET  formulation  is  given  by  Freeman  (1978). 

The  mean  flow  solution  described  above  will  respond  to  a  perturbation 
in  accordance  with  the  Navier-Stokes  equations  (1);  the  perturbation  may 
grow,  diminish  or  maintain  constant  amplitude  depending  on  both  the  nature 
of  the  mean  flow  and  of  the  perturbation  itself.  The  perturbation  is 
assumed  to  have  a  small  enough  amplitude  that  the  governing  equations  may 
be  linearized  (neglecting  terms  of  higher  order  than  the  first  in  the 
perturbation  amplitudes).  While  the  perturbation  itself  is  a  function  of 
horizontal  position,  as  well  as  altitude  and  time,  the  coefficients  of  the 
governing  equations  do  not  depend  on  horizontal  position.  For  this  reason, 
the  perturbation  may  be  assumed  to  be  a  sinusoidal  function  of  the  horizontal 
coordinate,  x. 
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-  =  .0  (z,_k,t)ei-‘-  3- (4) 

where  k_  is  the  vector  horizontal  wavenumber  of  the  perturbation;  a 

perturbation  of  more  general  horizontal  dependence  may  be  constructed  by 

superposition  of  these  solutions.  For  the  assessment  of  stability  the 

wavenumber  vector  k^  =  {k  ,  k  )  constitutes  an  additional  parameter  which 

x  y 

must  be  sampled. 

The  time  dependence  of  the  solution  v0(z,t)  determines  the  stability 
of  the  perturbation.  In  general,  the  governing  equations  form  an  initial 
value  problem;  starting  from  a  given  amplitude,  the  solution  will  change 
with  time  in  accordance  with  the  perturbation  equations  themselves.  The 
time  dependence  is  simplified,  however,  if  the  coefficients  of  the 
perturbation  equations  are  independent  of  time;  such  is  the  case  when  the 
solution  for  the  mean  flow  forms  a  steady  state,  as  in  the  case  of  the 
classical  Ekman  spiral.  In  general,  the  mean  flow  solution  is  not  steady 
(due  to  diurnal  oscillations,  inertial  oscillations,  initial  state  of  un¬ 
stable  heating,  etc.),  but  it  is  sufficiently  convenient  that  v/e  will 
make  this  assumption  for  the  interim.  Then  we  can  write. 


*0(z,k,t)  =  ^00(z,k)e-iu,t  B- ( 5 ) 

and  if  we  wish,  the  more  general  solution  may  be  constructed  by  super¬ 
position. 

We  now  carry  out  the  derivation  of  the  linearized  equations  govern¬ 
ing  the  perturbation.  By  subtraction  of  B-(2)  from  B-(l)  and  using  B- (4 ) 
we  obtain 
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,2  VW  ,  /  „H  \ 
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i  k/£  +  wz  =  0  .  B-(6) 

In  the  perturbation  equation  for  w  the  small  term  involving  f*  was  neg¬ 
lected  and  the  hydrostatic  approximation  B-(3)  was  used  to  eliminate  the 
vertical  component  of  the  mean  pressure  gradient. 

The  equations  B-(6),  together  with  boundary  conditions,  constitute 
the  linear  perturbation  problem  for  an  incompressible  fluid  in  a  very 
general  form.  The  coefficients  may  be  quite  arbitrary  functions  of 
altitude  and  time.  While  we  plan  to  return  to  this  general  problem  in  the 
future,  for  the  interim  we  assume  that  the  coefficients  do  not  depend  on 
time,  so  that  the  solution  form  B-{5)  is  applicable.  The  frequency  u 
constitutes  a  complex  eigenvalue  whose  value  is  determined  by  requiring 
the  solution  for  a  given  wavenumber  vector  k  to  satisfy  the  prescribed 
boundary  conditions.  For  each  value  of  there  is  a  multiplicity  of 
eigenvalues;  they  are  distinguished  from  each  other  by  a  mode  index 
corresponding  to  the  number  of  nodes  of  the  solution  in  the  altitude 
range. 

Let  us  now  substitute  the  time  dependence  of  B- ( 5 )  into  B-(6)  and 
rearrange  the  equations  for  easier  solution.  The  first  two  equations  of 
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B-(6)  are  combined  to  form  components  parallel  to  k_  (y  =  k/u)  and 

perpendicular  to  k  (n  s  k  u  -  k  v) 

y  x 

2 

(<  u2)z  -  W  y2  +  [i(u  -  k-U)  -  k2  kJJ]  y  -  k-TJ2  w  -  P  -  f'k  w 

P  X 

"  f  Cn  +  -  (ky  Ug  '  kx  V  =  0  » 

(Kv  Vz  "  “  n2  +  -  JL*D  -  1(2  #  n  -  (ky  U  -  kx  V)z  w  -  f'ky  w 

+  f  [y  +  z  k-U  p]  =  0  .  B-(7) 

p 

These  equations  are  coupled  together  only  by  the  Coriolis  term.  The 
equation  for  *n  does  not  contain  the  pressure. 

It  is  now  convenient  to  eliminate  the  pressure  from  the  w  equation 
by  using  the  first  of  B- ( 7) ;  we  also  eliminate  y  by  using  the  last  equation 
of  B- (6) .  When  we  do  this  the  order  of  the  resulting  equation  is  raised 
to  the  fourth.  Rather  than  raise  the  order  still  higher,  we  solve 

simultaneous  equations  for  n,  p  and  w.  The  equations  are  (neglecting  f ’ ) 
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B-(8) 


To  summarize,  we  have  derived  perturbation  equations  for  a  boundary  layer 
in  which  the  unperturbed  state  contains  an  arbitrarily  stratified  density 
and  wind  shear.  Large-scale  processes  are  represented  by  the  geostrophic 
wind  and  the  horizontal  divergence  of  the  mean  wind.  We  have  included 
Coriolis  terms,  but  have  assumed  incompressibility.  With  these  effects 
taken  into  account  the  equations  are  considerably  more  general  than  have 
been  studied  before.  They  include  as  special  cases  free  convection, 
internal  wave  propagation,  Ekman  spiral  instability,  etc.  If  the  Coriolis 
terms  are  omitted  the  n-equation  need  not  be  solved  simultaneously  with 
those  for  w  and  o.  The  retention  of  geostrophic  and  divergence  terms  does 
not  greatly  complicate  the  equations.  The  diffusion  terms  increase  the 
order  of  the  equations;  in  general,  the  equivalent  order  of  a  single 
equation  is  the  eighth,  but  if  the  Coriolis  term  is  ze^o  the  order  is 
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recuced  to  sixth.  The  constant  density  !neutral  stratification)  problem 
is  of  sixth  order,  and  becomes  of  fourth  order  when  f  =  0.  Omitting 
ciffusion  terms,  the  equations  are  of  fifth  order  unless  W  =  0  when  they 
are  of  second  order.  For  this  last  case  the  equations  are  readily  solved 
by  iteration  of  the  solution  of  the  second  order  equations  as  in  ZMODE 
(Milder,  1973).  It  is  thus  apparent  that  the  internal  wave  propagation 
problem  (which  is  of  second  order)  can  be  generalized  to  include  geo- 
strophic  terms  and  unstable  atmospheric  stratification,  as  well  as  the 
effect  of  wind  shear.  Numerical  solution  of  the  equations  could  be 
accomplished  by  several  different  methods,  of  which  the  iterative 
integration  method  used  in  ZMODE  has  been  alluded  to.  The  method  of 
time  integration  has  also  been  mentioned.  Each  of  these  methods  presents 
some  difficulties;  the  iterative  method  would  require  generalization  to 
higher  order  coupled  equations  as  would  the  eigenvalue  search  procedure, 
while  the  time  integration  method  suppresses  modes  with  smaller  eigen¬ 
values.  In  the  following  sections  a  matrix  solution  technique  is  outlined. 
This  method  has  the  advantage  that  all  eigenvalues  are  obtained  without 
iteration,  and  accuracy  is  readily  controlled. 
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B. 2  Difference  Equations 

The  matrix  method  of  solution  of  B- (8 )  requires  that  difference 
equations  on  a  finite  grid  of  points  replace  the  differential  equations. 
When  this  has  been  accomplished  the  resulting  system  of  linear  coupled 
algebraic  equations  can  be  expressed  as  a  matrix  equation  for  the  vector 
consisting  of  the  unknown  values  of  the  three  quantities  r,  o  and  w  at 
each  of  the  altitude  positions.  Since  the  frequency  u  enters  into  certain 
terms  of  these  equations  as  a  linear  factor,  the  complete  set  of  equations 
can  be  cast  in  the  form  of  a  generalized  eigenvalue  problem.  In  general, 
the  elements  of  the  coefficient  matrices  are  complex  and  the  resulting 
eigenvalues  may  be  complex.  The  real  part  of  the  eigenvalue  measures 
the  phase  speed  of  the  disturbance,  while  the  imaginary  part  measures  its 
growth  or  decay  rate.  If  the  atmosphere  is  divided  into  I  cells  or  grid 
points  the  values  of  the  unknown  eigenvector  X  (consisting  of  n,  p  and  w) 
are  31  in  number.  The  difference  equations  corresponding  to  B- (8)  are  of 
the  form 


AX  =  Bu>X 


B-(9) 


where  A  and  B  are  square  matrices  of  order  31.  The  exact  form  of  the  co¬ 
efficient  matrices  depends  on  the  difference  equations,  on  the  order  of 
arrangement  of  the  elements  of  the  eigenvector,  and  on  the  boundary 
conditions.  In  order  to  simplify  the  matrices  (and  facilitate  partitioning 
into  Ixl  submatrices)  the  order  of  the  eigenvector  is  chosen  to  be 

X  =  (Hi ,  n2  •  •  •  hj ,  c-j ,  ■  •  •  >.'j ,  Wj ,  W2  •  •  ■  Wj } 
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In  this  section  the  difference  approximation  to  the  differential  expressions 
will  be  developed.  We  assume,  at  this  time,  that  the  independent  variable 
X  is  divided  into  equal  intervals  IX  sucn  that  X  -  =  IAX;  the 

center  position  of  each  interval  is  associated  with  the  values  rj . ,  and 
w..  In  the  following  section  we  consider  the  transformation  from  the 
altitude  Z  to  X,  such  that  the  interval  in  Z  is  non-uniform,  as  required 
for  accuracy  in  a  planetary  boundary  layer  calculation.  Referring  to 
B-(8),  we  note  that  the  highest  order  of  derivative  is  tne  fourth,  which 
requires  a  minimum  of  five  adjacent  points  for  representation.  This 
implies  that  at  least  one  of  the  Ixl  submatrices  will  be  pentadiagonal 
(containing  non-zero  values  on  the  principal  diagonal  and  two  bordering 
diagonals  above  and  two  below  it).  We  elect  to  represent  lower  order 
derivatives  with  a  fourth  order  difference  expression  as  well,  since  to 
do  so  will  presumably  increase  the  accuracy,  will  invoke  a  consistent 
order  of  difference  approximation,  and  will  not  increase  the  number  of 
non-zero  matrix  elements.  Boundary  conditions  will  affect  the  difference 
expressions  at  points  immediately  adjacent  to  the  top  and  bottom  boundaries 
and  points  one  zone  removed  from  them.  In  the  interior  of  the  mesh  the 
fourth  order  centered  difference  approximations  to  the  first  four 
derivatives  are 
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B- (10) 


where  the  derivatives  are  evaluated  at  x.  and  the  subscript  on  the  dependent 
function  f.  denotes  the  position  at  which  it  is  evaluated.  The  above 
expressions  may  be  represented  more  concisely  in  matrix  notation 
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V.nile  it  is  quite  straightforward  to  derive  the  expressions  given  in 
B-(10),  it  is  desirable  to  forr.ulate  a  general  method  for  doing  so  which 
is  also  applicable  at  boundary  positions.  Consequently,  the  function  and 
its  derivatives  are  expressed  in  terms  of  neighboring  values  of  the 
function.  The  quartic  approxiration  is 


f(x)  =  f.  +  I  a.(x-x.)j  , 

j=l  J  1 


in  terms  of  which  the  derivatives  at  x^  are 


f<k>.K!ak  . 


We  can  now  find  the  four  unknown  coefficients  a,  by  evaluating  the  function  at 

J 

x  =  xi  +  nAx  (the  four  neighboring  positions). 


4  f^-  i  i 

fHn =  fi  +.E,  ir 4x0 "  • 


B- (12) 


where  n  =  -2,  -1 ,  0,  +  1 ,  +  2.  We  have  replaced  ai  by  the  corresponding 
cerivative,  and  now  have  five  equations  *or  the  function  and  its  four 
derivatives.  When  written  as  a  matrix  equation  B- ( 1 2 )  becomes 


The  inverse  of  C  gives  the  desired  coefficient  matrix  C"^  of  B-(ll). 

It  is  now  easy  to  write  the  C  matrix  by  inspection;  the  C-1  matrix 
is  obtained  by  numerical  inversion,  and  subsequent  operations  are  performed 
as  matrix  and  vector  products,  resulting  ultimately  in  the  elements  of  the 


t  A  and  B  matrices. 
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B.  3  Coordinate  Transformations 

The  finite  difference  intervals  of  the  altitude  Z  are  chosen  as  a 
function  of  altitude  in  the  SIGMET  bounaary  layer  computer  code.  Since 
the  mean  wind  shear  is  large  near  the  surface,  higher  resolution  near 
the  surface  results  in  greater  computational  efficienc...  An  analytic 
expression  has  been  chosen  which  contains  a  linear  and  a  logarithmic  term, 
such  that  near  the  surface  the  wind  velocity  obeys  the  law  of  the  wall. 

Our  application  of  this  transformation  is  based  on  the  form  chosen  for 
the  SIGMET  computer  code;  it  contains,  in  fact,  two  transformations.  The 
first  introduces  a  scaled  pressure  coordinate,  a,  the  purpose  of  which  is 
to  simplify  boundary  conditions  in  a  multidimensional  primitive  variable 
calculation;  for  the  present  application  this  transformation  is  not 
essential  and  introduces  some  complication  which  would  otherwise  be 
absent.  The  second  transformation  generates  a  non-uniform  mesh  in  the 
r  coordinate  in  order  to  improve  resolution  in  the  surface  boundary  layer; 
the  latter  transformation  is  useful  for  the  perturbation  equations  as  well 
as  the  equation  of  mean  flow.  To  maintain  compatibility  with  SIGMET  we 
retain  both  of  these  transformations  which  are  described  below. 

The  sigma  coordinate  is  defined  in  terms  of  the  mean  pressure  F, 


where  m  *  Fg  -  Fj  ,  and  Fg  =  mean  pressure  at  the  surface,  and  Pj  =  mean 
pressure  at  the  top  of  the  computational  region.  The  r  coordinate  has  the 
range  0  <  c  <  1,  where  zero  corresponds  to  the  top  of  the  mesh  and  unity 
to  the  surface.  Using  the  hydrostatic  relation  B-(3)  we  obtain 
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dc 

dz 


B- (14) 


The  transformation  to  a  non-uniform  c  grid  uses  a  log-linear  relation 


x  =  l-o  +  d  in  (l-o  +  7  ) 

o' 


where  a  and  oQ  are  parameters  determining  the  extent  and  severity  of  the 
non-uniformity.  We  will  require 


da  u 


l-o  i- 


B- (15) 


a£n[o  /(1+a  )]-l 

The  x-coordinate  will  be  divided  into  uniform  intervals  ^ - - - 

We  now  wish  to  express  the  Z-derivatives  (up  to  order  four)  which  appear 
in  B-(8)  in  terms  of  x-derivatives;  due  to  the  uniform  interval  on  X  we 
may  use  the  difference  expressions,  derived  in  Section  B-2,  to  evaluate 
the  X-derivatives. 


d  _  do  dx  d  _  gp  ,,  a  \  d  B- ( 1 6) 

dz  "  dz  da  dx  tt  k  1-a  +  oQ'  dx 

Taking  account  of  the  dependence  of  the  coefficient  of  B- ( 1 6 )  on  the 
higher  derivates  are  obtained  by  successive  application  of  the  above 
formulae. 
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Clearly,  equations  B- ( 1 6 )  and  B- ( 1 7 )  can  be  expressed  as  a  matrix,  relating 
the  z-derivati ves  with  the  x-deri vatives ,  having  the  form 
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where  the  coefficients  a.,  are  obtained  in  an  obvious  way  from  B-(16)  and 

J  ^ 

B-(17),  and  f^  5  — r  as  before.  Denoting  the  matrix  of  B-(18)  by  D,  the 


f ^ >  vector  can  be  evaluated  from  B-(ll)  to  give 
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Finally,  B- ( 1 9 )  can  be  multiplied  by  the  vector  consisting  of  the  co¬ 
efficients  of  the  derivatives  of  each  of  the  variables  q,  p,  w  in  the 
three  equations  of  B-(8);  the  results  are  the  vectors  giving  the  elements 
on  the  pentadiagonals  of  the  submatrices.  For  illustration,  consider  the 

simple  example  of  the  terms  from  the  third  equation  of  B-(8)  containing 

o 

n:  if  (q  +  —  n).  The  coefficients  of  q.  are  given  by  the  vector  result- 
p  p7  1 

from  (if  —  ,  if,  0,  0,  0}  DC  .  These  coefficients  are  a 

P 

function  of  altitude  due  to  the  height  dependence  of  pz_ 

The  above  procedure  is  used  to  obtain  all  of  the  matrix  coefficients 
in  the  interior  of  the  matrix.  Boundary  conditions,  however,  must  be  in¬ 
voked  on  zones  near  the  boundaries.  These  affect  the  coefficients  result¬ 
ing  from  the  two  zones  adjacent  to  the  top  boundary  and  the  two  adjacent 
to  the  bottom  boundary.  These  are  derived  in  the  next  section. 
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In  Section  2  difference  approximations  were  derived  between  the 
derivatives  entering  the  perturbation  equations  B- (8)  and  the  values  of 
the  function  at  mesh  points.  These  relations  involve  values  of  the 
function  two  intervals  above  and  below  as  well  as  at  the  cell  in  question. 
Clearly,  at  cells  adjacent  to  or  one  cell  removed  from  the  boundaries 
the  difference  expression  B-(ll)  is  not  suitable  because  it  would  require 
undefined  values  of  the  functions.  At  these  locations,  however,  the 
differential  equations  also  require  special  treatment  in  the  form  of 
boundary  conditions  describing  physical  conditions  of  the  interaction  of 
the  fluid  and  the  boundary.  These  boundary  conditions  are  often  idealized 
and  simplified,  corresponding  in  some  cases  to  the  necessity  of  truncating 
a  large  region  with  an  artificial  boundary  (such  as  the  upper  boundary 
of  the  planetary  boundary  layer). 

In  this  section  we  derive  relations  analogous  to  B-(ll)  which  are 
to  replace  B-(ll)  in  the  boundary  cells.  This  is  accomplished  by  modify¬ 
ing  the  equations  B- { 1 2 )  which  define  the  coefficients  of  the  quartic 
(and  which,  in  turn,  are  proportional  to  the  first  four  derivatives).  Each 
equation  corresponding  to  an  evaluation  of  the  quartic  outside  the  region 
is  replaced  by  an  equation  corresponding  to  a  boundary  condition.  It  is 
clear  from  this  construction  that  each  boundary  cell  will  require  two 
independent  boundary  relations,  while  cells  removed  from  the  boundary  will 
require  a  single  relation. 

The  boundaries  are  assumed  to  be  rigid.  At  the  surface  (the  lower 
boundary)  a  no-slip  condition  of  zero  velocity  is  invoked,  since  the  fluid 
is  assumed  to  be  viscous.  This  condition  also  applies  to  the  perturbation 
velocity,  and  implies  that  w  e  w  a  r.  a  0.  At  the  upper  boundary  we 


assume  the:  the  vertical  velocity  vanishes,  and  that  the-'e  is  no  shear 
stress  or  turbulent  transfer  across  it.  Applied  to  the  perturbation 
quantities  this  gives  w  =  wzz  =  =  -  =  (  We  also  assume  for  the 

interim  that  p  =  0  at  z  =  0.  This  set  of  boundary  conditions  is  sufficient 
for  the  differential  equations  (four  concitions  for  the  fourth  order  w 
equation,  and  two  conditions  for  each  of  the  second  order  -  and  p  equations). 
The  difference  equations,  however,  are  of  fourth  order  in  all  quantities; 
consequently,  some  non-physical  boundary  conditions  are  required  to  furnish 
the  additional  constraints  on  the  c  and  -  equations. 

Let  us  consider  the  w  equations,  where  we  have  w=wz=0atz=0 
and  w  =  w  =  0  at  z  =  H.  At  the  cells  adjacent  to  z  =  0  and  z  =  H  both 
of  the  boundary  conditions  are  applied  at  the  positions  of  the  boundary 
one-half  cell  thickness  away.  In  the  cells  one  removed  from  z  =  0  or 
z  =  H  only  one  of  the  boundary  conditions  is  needed;  we  choose  the 
condition  w  =  0  to  be  applied  at  the  bottom  and  wzz  =  0  at  the  top.  At 
the  lower  boundary  cell,  for  example,  the  equations  derived  from  the 
quartic  consist  of 

4  .  . 

f,  =  f,  +  I  a.  AxJ  2J  , 

O  1  1  J 

4  .  . 

f5  =  f,  +  Z  a.  AxJ  1J  , 

*  1  1  J 

4  . 

f0  -  0  =  f1  +  I  a.  Lx3  (-1/2)J  , 

f  =  0  =  I  ja.  Lx3'  (- 1 /2) 1  ,  B-(20) 

1  J 

where  the  first  two  equations  correspono  to  evaluating  tne  functions  in 
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the  mesh  interior  at  2 lx  and  lx  above  the  center  of  the  boundary  cell, 
the  third  equation  corresponds  to  w  =  0  at  -1/2  lx,  and  the  fourth 
corresponds  to  wz  =  0  at  -1/2  lx.  Equations  B-(20)  can  be  written  in  the 
matrix  form  B- ( 1 3 ) 

(Bottom  Boundary) 
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B-(21) 

The  remaining  boundary  matrices  are  obtained  in  a  similar  way.  Suppress¬ 
ing  the  second  matrix  factor  of  C,  which  is  the  same  for  each  case,  we 
have 


(Next-to-Bottom  Boundary) 
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B- ( 24 ) 

In  the  calculations  which  have  been  carried  out  to  date  we  have  used  these 
same  matrices  for  the  boundaries  of  the  --and  o-submatrices ,  as  well  as  for 
the  w-submatrices.  It  will  be  desirable  to  reexamine  the  z  boundary 

conditions  and  to  relax  the  condition  of  :  =  0  at  the  top  boundary.  In 

order  to  do  so  it  will  be  necessary  to  define  separate  C-matrices  at  the 
boundaries  for  n,  o  and  w.  It  is  likely  that  the  boundary  conditions  at 

the  top  of  the  mesh  will  not  turn  out  to  be  very  significant  because  the 

perturbation  amplitudes  should  become  small  at  high  altitude.  If  the 
boundary  is  located  high  enough  the  solutions  should  not  be  affected  by 
conditions  there. 


» 


» 


5.5  Numerical  Solution  of  Matrix  Eigenvalue  Equat'c 

We  employ  a  standard  subroutine  package  for  solving  the  complex 
eigenvalue  problem.  In  order  to  place  B-(9)  in  t^e  ^orm  required  by  the 
solver  we  multiply  by  the  inverse  of  the  3  matrix 


B_1AX  =  A 


B- (25) 


The  B  matrix  inverse  was  obtained  quite  economically  since  it  has  a  block 
diagonal  form 


o 


B  = 


0  B22  0 


V 


0  0  I 


,  where  I  is  the  identity  matrix 


/ 


and  the  element  B22  is  purely  imaginary.  We  form  and  construct 


22 
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,-i 


^-i 


22 
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1 


■1 


The  matrix  B  A  is  then  submitted  to  the  eigenvalue-eigenvector  subroutine. 

The  computer  solution  of  the  complex  general  eioensystem  was  per¬ 
formed  with  routines  from  the  EISPACK  (Smith,  1976)  mathematics  library.  These  are 
FORTRAN-based  routines  published  under  the  auspices  of  the  Applied 
Mathematics  Division,  Argonne  National  Laboratory. 

The  eigensystem  package  includes  a  driver  routine  "CG"  which  calls 
a  recommended  sequence  of  subroutines  to  find  the  eigenvalues  and  (if  re- 
auested)  eigenvectors  of  a  complex  general  matrix. 
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The  sequence  of  routines  called 


(a)  CBAL:  balances  a  complex  general  matrix. 

(b)  CORTH:  given  the  balanced  complex  general 

matrix,  reduces  it  to  a  sub-matrix 
of  upper  Hessenberg  form  by  unitary 
similarity  transformations. 

(c)  COMQR:  (eigenvalues  only)  or  C0MQR2  (eigen¬ 

values  and  eigenvectors):  finds  the 
eigenvalues  and  eigenvectors  of  the 
complex  upper  Hessenberg  matrix  using 
the  QR  algorithm  of  Francis  (1961,  1962). 

(d)  CBABK2:  forms  the  eigenvectors  of  the  original 

complex  general  matrix  by  back- 
transforming  those  of  the  corresponding 
balanced  matrix. 

For  convenience  in  comparing  with  test  results  the  eigenvalues  may 
be  presented  in  two  alternative  formats;  as  real  and  imaginary  parts  of 
the  eigenvector,  or  as  the  magnitude  (amplitude)  and  phase  of  the  eigen¬ 
vector.  Since  each  eigenvector  contains  an  arbitrary  complex  amplitude, 
we  have  normalized  the  magnitude  to  1  for  presentation,  and  require  the 
phase  <i>  to  be  in  the  range  -  v  £  <f>  <_  m.  Examples  of  this  presentation  are 
given  in  subsequent  sections. 

Based  on  comparisons  discussed  in  Section  B-6,  the  computer  code 
and  eigenvalue  subroutine  are  found  to  give  quite  accurate  estimates  of 
the  lower  mode  eigenvalues.  However,  we  have  found  that  the  31  x  31  ele¬ 
ment  eigenvalue  subroutine  fails  under  a  limited  range  of  conditions.  We 
believe  that  these  conditions  correspond  to  double  roots  of  the  indicia! 
equation,  which  can  arise  under  certain  idealized  conditions.  This 
problem  was  encountered,  for  example,  when  f  =  c  =  0.  We  have  shown  that 
under  these  circumstances  the  indicial  ecuetion  becomes 

(A-j  1  -  oil)  (B'2\  A22  -  tjl)  (A33  -  .  I)  =  0 


.-.'here  each  of  the  quantities  is  an  I  x  I  matrix.  For  tre  cited  case. 

A^  =  A^,  and  each  of  I  eigenvalues  is  ^eoeated.  In  tris  instarce  it 
is  possible  to  reduce  the  order  of  the  equations,  for  exar.pl e  by  '"emo.ing 
the  equations  forr,  which  are  decoupled  from  the  others.  The  resulting 
21  x  21  element  problem  is  solved  by  the  standard  technique  without 
difficulty.  When  parameters  are  chosen  to  correspond  to  a  general  problem 
the  difficulty  of  eigenvalue  degeneracy  is  not  encountered.  Consequently, 
it  appears  that  the  failure  of  the  technique  will  be  restricted  to  a 
limited  range  of  parameters  where  a  coefficient  (such  as  f  or  o’  )  is  very 
small ,  but  not  zero. 

We  anticipate  that  difficulty  also  right  occur  associated  with  a 
critical  level  (where  u  =  K*U)  but  have  not  examined  this  point  in 
detail  (see  Section  D).  It  is  expected  that  the  amplitude  of  the  perturb¬ 
ation  will  increase  near  the  critical  level,  causing  a  breakdown  of  the 
equations.  However,  viscous  damping  will  act  to  limit  the  growth;  it 
is  likely  that  a  careful  inclusion  of  absorption  and  choice  of  zoning  will 
eliminate  difficulties  associated  with  critical  layers. 


Test  Calculations 


r 


* 


The  computer  code  described  in  lections  B-l — 5-5  above  res  teen 
programmed  and  tested  against  numerical  solutions  obtained  with  ether 
techniques  and  by  other  investigators.  These  comparisons  verify  that  the 
coding  is  substantially  correct  and  tnat  solutions  are  being  obtained  with 
reasonable  efficiency.  Two  classes  of  test  problems  were  evaluated  and  com¬ 
parisons  were  made. 

In  the  first,  the  solutions  of  the  inviscid  negative  buoyancy  stability 
problem  were  obtained.  These  can  be  corpared  with  eigenvalues  calculated 
with  the  ZMODE  iterative  integration  corputer  code,  which  we  have  employed 
previously.  When  there  is  no  wind  these  solutions  correspond  to  internal  wave 
propagation  with  real  eigenfrequencies.  Using  an  analytic  profile  of  buoyancy 
and  zero  wind  speed  we  calculate  the  eigenvalues  as  a  function  of  wave  number 
using  different  numbers  of  zones  in  the  vertical  direction.  We  found  that 
the  eigenvalues  were  all  real,  corresponding  to  undamped  propagating  internal 
waves  and  that  they  occur  in  positive-negative  pairs,  corresponding  to  propaga¬ 
tion  in  a  given  wavenumber  direction  anc  in  a  direction  180°from  the  first. 

In  these  comparisons  uniform  zoning  was  used.  For  a  case  in  which  the  buoyancy 
frequency  profile  is  given  by 


> 


N2  =  C1  X  4exp(-2X) ,  0  -  Z  -75., 
where  Z  =  is  the  vertical  height, 


> 


X  =  -Z/C3  and  Cj  =  3.41  X  10*",  5-  =  25.0. 


» 
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.-.e  have  compared  both  eigenvalues  and  eigenfunctions,  “his  case  served  as  a 
verification  that  the  matrix  eigenvalue  subroutines  that  we  employ,  after  some 
initial  difficulty,  are  calculating  correct  eigenfunctions.  Using  20  zones 
we  obtained  the  following  eigenvalues  for‘  the  case  ;k'  =  2.5cyc/km,  (we 
compare  them  with  the  corresponding  values  resulting  from  the  ZMODE  computer 
code ) : 


Mode 

PERT 

ZMODE 

1 

.00310 

.00311 

2 

.00145 

.00146 

3 

.00095 

.00095 

4 

.00070 

.00074 

We  also  display  the  eigenfunctions  corresponding  to  the  vertical 
velocity  perturbation  w,  as  calculated  by  the  two  methods.  In  Fig.  B-l 
the  eigenfunctions  for  mode  1  are  shown;  it  was  necessary  to  rotate  the  eigen¬ 
function  by  an  arbitrary  phase  and  normalize  it  to  unit  maximum  amplitude. 

Figs.  B-2  and  B-3  correspond  to  the  second  and  third  mode  numbers.  In  each 
of  these  cases  the  phase  of  the  eigenfunction  is  independent  of  altitude. 

These  comparisons  indicate  that  excellent  agreement  between  PERT  and  ZMODE 
eigenvalues  and  eigenfunctions  is  achieved.  Agreements  for  higher  modes 
oecomes  less  quantitative  as  the  mode  number  approaches  the  number  of  vertical 
zones.  It  should  be  feasible,  however,  to  form  arbitrarily  large  mode  number 
solutions  by  increasing  the  number  of  zones. 

The  second  test  problem  corresponds  to  the  instability  of  Ekman 
ooundary  layer  flow.  An  idealized  version  of  this  problem  was  studied  by 
_illy  (1965),  who  assumed  that  the  atmosohere  is  neutrally  stable  and  that  the 
cirfusivities  are  independent  of  altitude.  The  wind  field  is  a  solution  of  the 
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►  steady  flow  equations  for  a  time-and  altitude-independent  geostrophic  wind, 
and  corresponds  to  the  classical  Ekman  st'-'al. 

►  U  =  V  (l-e‘Z/D  cos  2/D) , 

V  =  vn  e”Z/D  sin  Z/D, 

b 

I 

« 

where  VQ  is  the  height-independent  geostrophic  wind,  and  D  =  (v/H)1  is  the 
Ekman  depth.  Lilly  calculated  stability  diagrams  for  the  Orr-Sommerfeld 

*  equation  and  for  a  more  general  problem  in  which  the  Coriolis  terms  are 
included  in  the  perturbation  equations.  Our  test  consisted  in  comparing  with 
a  large  number  of  cases  reported  by  Lilly.  For  the  Orr-Sommerfeld  problem, 

*  parameters  close  to  the  point  of  maximum  instability  were  selected  for  several 

V  D 

values  of  the  Reynolds  number  R  =  where  is  the  geostrophic  wind  speed, 
v  is  the  altitude-independent  turbulent  diffusivity,  and  D  is  the  characteristic 

*  Ekman  length  (see  Lilly,  1966  for  a  comprehensive  discussion  of  this  problem). 

In  constrast  to  the  internal  wave  problem,  the  eigenfunctions  of  the  Ekman 
problem  depend  strongly  on  the  direction,  as  well  as  the  magnitude,  of  the 

®  wavenumber  of  the  perturbation.  We  examined  several  cases  in  which  the 

direction  of  the  wavenumber  vector,  the  Reynolds  number,  and  the  zoning  were 
varied.  In  agreement  with  Lilly's  results,  we  found  that  R  =  65  is  stable 

*  and  that  R  =  110  and  R  =  500  are  unstable.  In  both  of  the  latter  cases, 
instability  is  due  to  a  single  mode  of  the  perturbation;  all  other  modes  are 
strongly  damped.  We  also  compared  the  real  and  imaginary  parts  of  the  eigen- 
frequency  with  Lilly's  calculations  (by  inspection  of  his  small  diagrams). 

We  found  general  agreement  with  the  growth  rate  and  phase  speed  of  the  dis- 

» 
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zurbarce  in  the  vicinity  of  maximum  instability  as  well  as  at  positions  where 
the  growth  rate  is  smaller  or  negative.  With  few  vertical  zones  (corresponding 
to  the  initial  calculations)  the  eigenvalues  did  not  compare  as  well  as  in  the 
internal  wave  calculations.  Increasing  the  number  of  uniform  size  zones 
improved  agreement,  but  substantial  differences  ('-  25  in  worst  cases)  still 
remained  with  32  zones  in  the  vertical.  Kodifina  boundary  conditions  in  zones 
one  removed  from  the  boundary  (to  relax  the  boundary  constraint)  restored 
agreement.  Since  the  wind  profile  contains  large  shear  near  the  surface,  a 
nonuniformly  zoned  mesh  was  also  tested.  ~his  produces  some  improvement,  but 
we  have  not  had  the  time  to  determine  optimum  mesh  parameters  or  the  enhance¬ 
ment  of  accuracy.  Subsequently,  the  Coriolis  terms  were  added  to  the  matrix 
equations  and  a  test  of  the  complete  system  of  equations  was  made.  We  calculated 
several  cases  recorded  by  Lilly  corresponding  to  the  wavenumbers  where  the 
inviscid  and  the  parallel  instabilities  have  maximum  values.  In  addition  to 
the  eigenvalues  we  were  able  to  compare  our  eigenfunctions  with  those  given 
by  Lilly.  For  a  Reynolds  number  of  R  =  110  the  local  maximum  value  of  in¬ 
stability  due  to  the  parallel  mode  is  located  near  k  =  0.3  and  e  =  -11°. 

For  these  values  of  the  parameters  we  obtain  =  .124,  wj  =  .0052,  which 
are  in  good  agreement  with  Lilly.  In  order  to  compare  eigenfunctions  we  formed 
the  absolute  value  of  the  normalized  amclituae  and  the  phase  of  the  eigen¬ 
functions.  In  Fig.  B-4  we  show  these  quantities  for  the  w  eigenfunction; 
they  are  in  good  agreement  with  Lilly.  The  inviscid  mode  has  its  max¬ 
imum  near  k 1  =  0.5,  e  =  8°;  the  corresponding  eigenvalues  are 
, p  =  .0461,  =  .0021,  and  the  eigenfurctions  are  shown  in  Fig  B-5.  As 

indicated  by  Lilly,  the  Coriolis  terms  -ocify  the  stability  of  the  boundary 
layer  significantly  through  the  paralle'  instablity  mechanism.  The  effect 


*  2*  these  terms  is  to  lower  the  critical  Re.-nclds  number-  and  to  cr.anae  the 
character  of  the  unstable  motion.  Consecuently ,  it  is  imoortant  to  include 
tr.ese  terms  in  the  stability  analyses  of  a  neutrally  stable  bouncary  layer 

t  even  though  at  first  consideration  one  would  anticipate  that  the  Coriolis 

terms  are  too  small  to  affect  the  results. 

It  should  be  remarked  in  passing  that  the  above  formulation  does  not 

>  include  all  of  the  terms  containing  the  Coriolis  parameter.  Referring  to  (8) 

we  find  that  the  perturbation  equations  for  *■  and  w  have  terms  proportional 

to  (U  ,  V  )  and  f  which  we  discuss  below, 
y  y 

*  On  the  basis  of  the  above  comparisons  it  is  now  quite  likely  that  the 
major  terms  of  the  perturbation  computer  codes  are  correct,  and  that  the 
numerical  methods  developed  for  its  solution  are  working  satisfactorily. 

►  There  remain,  however,  a  large  number  of  terms  which  have  not  been  tested. 

These  include  terms  involving  gradients  of  diffusivity,  vertical  mean  velocity, 
and  geostrophic  terms  in  the  perturbation  equations.  It  will  be  necessary 

*  to  test  these  terms  as  well  in  the  course  of  performing  a  comprehensive  study 
of  boundary  layer  instability. 

A  difficulty  arises  in  trying  to  explore  the  many  parameters  which  the 

*  computer  code  contains.  Ultimately  these  free  parameters  will  be  constrained 
by  the  realizability  of  the  mean  flow  solution.  That  is,  it  is  not  practical 
or  realistic  to  independently  vary  the  many  degrees  of  freedon  presented  by  the 

*  mean  flow  quantities,  since  they  must  be  related  to  each  other  in  such  a  way 
that  they  represent  a  consistent  solution  of  the  mean  boundary  layer  equations 

fas,  for  example,  in  SIGMET).  This  consideration,  in  fact,  constrained  our 

^  formulation,  dictating  several  features  of  PERT,  such  as  the  matrix  form  and 

toe  non-uniform  zoning.  Several  calculations  have  been  performed  in  this  mode, 
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and  they  are  reported  in  Section  B-7. 
carried  out  some  calculations  designed 
obtaining  new  boundary  layer  stability 
reported  in  the  next  section. 


however,  chronologically 
to  show  the  capabilities 
results.  These  findings 


we  first 
of  PERT  by 
are  also 
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I  B-7  Stability  of  Boundary  Layer  Flo;-. 

In  the  previous  section  tests  o'  the  PERT  computer  code  were  described 
for  two  classes  of  problems  differing  greatly  in  physical  content;  wave 
I  propagation  in  a  stably  stratified  medium,  and  secondary  circulation  of  an 

idealized  shear  flow  in  a  neutrally  stratified  atmosphere.  These  problems 
show  that  the  technique  gives  expected  results  in  these  cases,  so  that  we  may 
l  now  look  at  several  more  general  boundary  layer  problems.  In  this  section  the 

results  of  calculations  sampling  some  of  the  additional  parameters  which  affect 
boundary  layer  stability  are  reported.  The  combined  effects  of  shear  and 
I  stratification  appear  not  to  have  been  investigated.  cirst  we  explore  the 

effect  of  a  constant  density  gradient  added  to  the  Ekman  flow  problem  invest¬ 
igated  in  Section  B-6.  The  mean  wind  configuration  of  the  Lilly  study  is  retained, 
»  the  diffusivities  are  equal  and  independent  of  altitude  as  before,  and  the 

density  gradient  is  taken  to  be  independent  of  altitude  as  well.  We  will 
explore  how  the  Orr-Sommerfeld  stability  diagrams  are  modified  by  this  constant 
§  density  gradient;  we  expect  to  find  increased  instability  when  the  atmosphere 

is  positively  buoyant,  and  stabilization  of  the  shear  instability  for  sufficiently 
large  stable  stratification.  The  density  gradient  introduced  into  the  perturba- 

•  tion  equations  is  a  dimensionless  quantity  scaled  by  the  Ekman  layer  depth, 

D,  employed  in  the  neutrally  stable  case.  The  figures  and  discussion  below 

7T  r 
* 

refer  to  this  scaled  quantity  o2*  s  — - .  A  typical  value  of  the  potential 

^  -5  -1 

*  density  gradient  in  the  atmosphere  of  10  m  gives  a  scaled  density 
gradient  -  3  x  10’^. 

For  the  case  of  the  Orr-Sommerfeld  equation  with  R  =  110  we  have 

♦  explored  the  dependence  of  the  growth  rate  on  the  density  gradient,  oz*. 

Using  K  =  0,  K  =  0.5  (near  the  point  cf  maximum  instability  of  the  neutrally 
x  y 

stratified  case),  we  varied  the  density  gradient  as  shown  in  Fig.  B-6. 
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A  very  small  stable  density  gradient  is  sufficient  to  acnieve  stability 
(:  *  -  -  .0002);  more  strongly  stable  atmospheres  exhicit  damping  for  all 
modes  of  the  perturbation.  It  is  interesting,  however,  that  the  growth  rate 
does  not  continue  to  decrease  in  the  stable  region  because  another  mode  assumes 
the  role  of  the  largest  imaginary  eigenvalue.  This  moce  becomes  dominant  at 
c2*  =  -  .0004,  and  its  magnitude  decreases  slowly  for  more  stable  stratification. 
This  mode  (and  other  modes  having  about  the  same  damping  rate)  is  probably 
an  internal  wave  mode,  which  becomes  more  prominant  as  the  Brunt-Vaisala 
frequency  is  increased.  The  dependence  of  the  growth  rate  on  the  magnitude 
and  direction  of  the  wavenumber  at  oz*  -  -  .0001  for  R  =  110  is  shown  in  Fig. 

B-7.  This  shows  that  the  configuration  of  the  unstable  region  is  not  substantially 
changed,  but,  as  expected,  the  extent  of  the  unstable  region  is  decreased. 

The  character  of  the  solution  also  changes  rapidly  when  the  density 
gradient  becomes  positive.  In  Fig.  B-6  another  mode  becomes  most  unstable  when 
o2*  -  .0002.  For  this  mode  the  growth  rate  is  a  rapidly  increasing  function  of 
P2*.  For  typical  values  of  unstable  atmospheric  stratification  the  predominant 
mode  will  be  this  convectively  unstable  one;  the  shear  instability  will  not 
influence  the  solution  appreciably.  It  appears  from  a  couple  of  calculations 
that  the  convective  mode  is  very  much  less  directional  than  the  shear  mode. 

We  have  also  performed  calculations  for  the  Orr-Sommerfeld  problem  with 

R  =  500,  for  which  the  unstable  region  in  K  ,  K  space  is  larger.  In  this 

X  y 

case,  as  indicated  in  Fig.  B-8,  the  density  gradient  must  be  more  negative 
(c2*  »  -  .0015)  to  overcome  the  shear  instability.  As  in  the  lower  Reynolds 
number  case,  there  is  a  crossing  of  mode  trajectories.  The  damping  rate  for 
the  most  nearly  unstable  mode  in  the  larger  gradient  region  is  now  quite  small, 
corresponding  to  the  larger  Brunt-Vaisala  •‘reouency. 


In  view  of  the  observed  effect  c£  tne  Coriclis  *orce  on  tre  Ekman 
stability  criterion,  we  have  repeated  several  of  the  ca'culations  assessinc 
the  effects  of  stratification  with  the  Coriolis  terms  included.  First,  we 
examine  the  effect  of  the  same  terms  retained  by  Lilly  and  included  in  our 
Coriolis  test  calculations  (see  Section  E-6).  The  effect  of  density  gradient 
on  growth  rate  for  this  case  is  shown  in  Fig.  B-9.  As  expected  from  Lilly's 
diagrams  the  instability  near  =  0  is  decreased  for  •  =  0.5,  =  0. 

However,  for  stable  stratification  the  "gravity  wave"  mode,  which  soon  dominates 
the  stability  diagram,  shows  very  little  dependence  on  the  Coriolis  parameter. 
Consequently  for  pz*  larger  than  that  for  mode  crossover,  the  most  unstable 
mode  growth  rate  approaches  that  with  f  =  0.  Similarly,  for  unstable  strati¬ 
fication  the  "convective"  modes  are  not  dependent  on  f  so  that  sufficient  in¬ 
stability  also  reduces  the  influence  of  the  Coriolis  paramiter.  These  properties 
are  displayed  in  Fig.  B-9,  by  comparison  with  Fig.  B-6  (with  f  =  0).  It  is 
interesting  that  the  mode  responsible  for  instability  at  :  z*  =  0  shows  markedly 
smaller  growth  rate  with  the  Coriolis  term  present  over  the  full  range  of  density 
gradient. 

The  effect  of  the  large-scale  convergence/divergence  of  the  mean  flow 
field  also  modifies  atmospheric  stability.  We  have  explored  this  effect  by 
adding  the  terms  in  W  and  D  to  the  neutral  Ekman  boundary  layer  perturbation 
equations.  Starting  from  the  Lilly  formulation  in  which  the  standard  Coriolis 
terms  are  included,  we  added  a  height-independent  divergence  term  D  (which 
may  be  either  positive  or  negative).  The  corresponding  mean  vertical  velocity 
is  1  inear  with  height  z, 

W  =  -  Dz. 
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Start’r.g  frcm  tr e  case  R  =  ill,  \  =  0,  =  0.5,  we  have  varied  0  tc  assess 

A  _/ 

’ ts  effect  on  the  growth  rate  of  the  most  unstable  eigenvalue.  In  Fic.  B- 10 
the  growtn  rate  is  a  function  of  the  scaled  divergence  in  a  range  of  values 
actually  found  in  the  atmosphere .  We  find  that  the  growth  rate  is  quite 
insensitive  throughout  this  range  of  convergence  and  divergence.  It  should 
be  remarked  that  the  effect  of  moisture  condenation  is  not  taken  into  account 
in  this  investigation;  the  effect  of  lifting  condensation  on  convective  instability 
will  be  much  larger. 

Finally,  we  have  performed  several  calculations  using  the  capability  to 
transfer  data  from  the  SIGMET  code  to  the  PERT  code.  In  this  procedure  all 
of  the  mean  variable  profiles  needed  to  initiate  the  linear  perturbation 
calculation  result  from  the  unperturbed  boundary  layer  calculation.  Data  are 
extracted  from  a  particular  time  of  a  time-dependent  calculation,  even  though 
this  is  not  entirely  consistant  with  the  time-independent  perturbation  for¬ 
mulation.  However,  the  SIGMET  solution  over  a  water  surface  is  slowly  changing. 
These  calculations  differ  primarily  from  previous  calculations  in  having 
profiles  of  mean  variables  which  are  not  simplified  for  the  sake  of  convenience. 

The  wind  field  can  change  in  nagnitude  and  direction,  the  buoyancy  field  may 
contain  simultaneously  megative  and  positive  regions,  and  the  turbulent  transfer 
coefficients  may  vary  with  altitude  as  required  by  the  boundary  layer  equations. 

Several  mathematical  difficulties  prevented  us  from  carrying  out  these 
calculations  until  very  date  in  the  investigation.  Consequently,  we  have  not 
been  able  to  explore  the  many  parameters  governing  this  boundary  layer  stability 
problem.  We  are  restricted  to  sampling  the  stability  diagrams  for  a  single  case. 

The  SIGMET  problem  has  oeen  reported  in  Section  P;  it  corresponds  tc  cycle 
12  (time  =  1608  local  time)  of  a  marine  summer  boundary  layer  having  a 


geostrophic  wind  speed  of  10  n  s.  Some  cf  the  profi’e  :eta  are  snown  in  Figs. 
A-l  and  A-2.  We  repeated  the  calculation  with  lower  retention  using  20  zones 
in  order  to  reduce  the  computational  expense  of  the  calculations.  Sub¬ 
sequently,  a  survey  was  made  of  linear  stability  as  a  *Jnct"'on  o*  wavenumber, 
"he  dispersion  diagram  for  these  calculations  is  shown  -in  Fig.  B-ll.  A 
number  of  the  modes  appear  to  be  internal  waves  displaced  by  the  geostrophic 
wind  speed.  The  up-wind  traveling  waves  are  more  complicated,  displaying 
several  mode  crossings.  These  later  waves  will  encounter  critical  levels  at 
lower  altitudes  of  the  boundary  layer.  Apparently,  a  shear-unstable  mode 
is  also  present.  The  magnitude  of  the  growth  rate  of  this  node  increases  with 
wavenumber  as  shown  in  Fig.  B-12.  The  similarity  of  the  dispersion  diagram 
with  that  from  the  internal  wave  study  of  Section  D  is  quite  striking.  Further 
investigation  of  this  an  similar  problems  is  warrented. 
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E. 8  Summary  and  Remarks 

The  linear  stability  of  boundary  layer  flow  has  been  investigated 
u°der  quite  general  conditions  of  wind  shear,  stratification  and  turbulent 
excnange.  Perturbation  equations  have  been  derived  admitting  arbitrary 
dependence  of  mean  flow  quantities  on  altitude,  and  these  are  formulated  for 
numerical  solution  as  a  matrix  eigenvalue  problem.  A  computer  program  (PERT) 
containing  this  formulation  have  been  developed  to  evaluate  quite  general 
boundary  layer  stability  problems.  Calculations  have  been  performed  dupli¬ 
cating  known  results  for  internal  wave  propagation  and  Ekman  flow  instability. 
The  stability  of  generalized  Ekman  flow  in  stratified  atmospheres  has  also 
been  studied.  We  have  determined  the  amount  of  negative  buoyancy  needed  to 
stabilize  the  shear  instability  and  have  displayed  some  of  the  systematic 
features  of  the  stability  diagram.  An  advantage  of  the  matrix  method  is  that 
all  of  the  lowest  modes  of  the  eigenvalue  problem  are  obtained.  Using  the 
resulting  data  it  then  becomes  possible  to  trace  the  competition  of  several 
modes.  For  the  problem  of  a  stratified  shear  boundary  layer  we  found  that 
the  crossing  of  different  mode  trajectories  is  responsible  for  the  changes 
in  stability  character  when  stratification  departs  from  neutral  in  both  the 
positive  and  the  negative  direction.  In  the  case  studied  we  found  that  only 
the  node  responsible  for  instability  near  neutral  stratification  is  sensitive 
to  the  Coriolis  term. 

We  have  also  investigated  the  effect  of  the  large-scale  convergence 
or  divergence  of  the  horizontal  wind  field  on  stability.  This  field  is  almost 
always  present  due  to  synoptic  weather  features  such  as  cyclones  and  anti¬ 
cyclones.  When  this  term  and  the  resulting  mean  vertical  velocity  were 
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added  to  the  neutrally  stable  Ekman  boundary  layer  we  found  that  the  growth 
rate  of  instability  was  modified,  but  that  within  the  range  of  realistic 
values  of  convergence  or  divergence  the  magnitude  was  not  appreciably  changed. 
We  conclude  that  its  effect  on  secondary  circulations  not  involving  moist 
processes  (which  have  not  been  included)  is  minimal. 

A  version  of  the  PERT  computer  code  has  been  developed  which  accepts 
data  from  the  one-dimensional  boundary  layer  computer  code,  SIGMET.  To  demon 
strate  this  capability  a  linear  stability  analysis  of  a  case  involving  a 
marine  atmosphere  has  been  performed.  The  selected  case  contains  stratifi¬ 
cation  which  is  slightly  stable  through  the  mixing  layer  and  quite  stable 
above.  Rather  small  turbulent  intensity  is  present  in  the  lower  layer. 

Several  modes  may  be  identified  as  corresponding  to  damped  internal  waves. 

One  mode  of  the  eigenvalue  spectrum  displays  exponential  growth;  this  pre¬ 
sumably  corresponds  to  the  shear  layer  instability  (without  Coriolis  effect, 
since  this  term  is  currently  omitted  from  the  version  of  PERT  used  in  this 
calculation) . 

In  Section  D  some  additional  PERT  calculations  are  presented  in 
connection  with  an  investigation  of  the  internal  wave  spectrum  to  be  expected 
in  a  marine  atmosphere. 

To  date,  we  have  been  able  to  indicate  qualitatively  the  range  of 
capabilities  of  the  linear  stability  analysis  through  the  illustrative  cal¬ 
culations  described  above.  Due  to  limitations  of  funds  and  time  we  have  not 
been  able  to  exploit  this  tool  to  its  full  capability.  In  assessing  our 
progress  toward  understanding  instability  and  secondary  circulations  of 
marine  atmospheres  it  is  important  to  place  linear  analysis  in  perspective. 


The  advantages  of  this  approach  are  quite  significant.  As  we  have  shewn, 
essentially  all  physical  effects  can  be  taken  into  account,  and  me  resulting 
calculations  make  quite  modest  demands  in  terms  of  comDuter  time.  Kitr'  this 
method  it  is  feasible  to  carry  out  hundreds  of  separate  calculations,  tnereby 
surveying  the  stability  parameter  space  more  thoroughly  than  possible  by  other 
methods.  The  accuracy  of  the  linear  perturbation  approach  is  more  difficult 
to  evaluate,  since  it  depends  on  the  particular  response  of  the  system.  To 
illustrate,  the  internal  wave  solutions  are  either  undamped  or  slightly 
damped  by  viscosity.  They  will  frequently  have  small  enough  amplitude  to 
satisfy  the  linear  approximation.  Under  certain  conditions,  however,  they 
will  grow  in  amplitude  until  they  can  no  longer  be  considered  to  be  of  small 
amplitude;  they  then  modify  the  "unperturbed"  solution  and  are  no  longer 
governed  accurately  by  the  linear  equations.  When  there  exists  an  unstable 
mode  we  expect  that  the  linear  solution  ultimately  will  break  down.  In  that 
event  the  linear  equations  suggest  the  conditions  for  the  onset  of  instability 
but  give  no  information  on  the  asymptotic  state  (if  one  exists)  of  the 
secondary  circulation.  Clearly,  there  is  no  substitute  for  a  nonlinear 
calculation  if  one  requires  information  about  the  state  that  the  system 
reaches  at  late  time.  Information  about  linear  and  nonlinear  behavior  appears 
to  be  complementary  in  many  respects. 

Our  linear  perturbation  studies  should  be  augmented  with  nonlinear 
calculations,  as  proposed  earlier.  But  it  will  also  be  desirable  to  extend 
the  linear  calculations  in  several  respects.  We  are  currently  improving  the 
accuracy  of  the  PERT  code  by  introducing  double  precision  arithmetic  in 
selected  parts  of  the  code.  There  is  also  potential  improvement  for  more 
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careful  boundary  conditions  and  zoning.  From  the  physical  point  of  view 
further  work  is  needed  on  time-dependent  unperturbed  solutions,  on  compressi¬ 
bility,  and  on  moisture  effects.  Additional  work  on  coupling  with  the  boundary 
layer  codes  will  enhance  our  ability  to  investigate  realistic  problems  in 
convective  instabil ity.  Finally,  as  mentioned  earlier,  we  believe  that  the 
applications  of  the  PERT  code  are  far  from  exploited.  Considerable  informa¬ 
tion  can  be  obtained  about  the  systematics  of  stability/instability  of 
marine  boundary  layers.  Not  the  least  of  the  benefit  will  come  from  enhanced 
understanding  of  the  dynamics  of  competing  modes. 
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NORMALIZED  EIGENVECTOR 

GURE  B-l.  Comparison  of  Mode  1  eigenfunctions  for  PERT  and  ZMODE  internal 
wave  calculations,  jLj  =  2.5  cyc/km. 
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FIGURE  B-7.  Growth  rate  of  most  unstable  mode  as  a  function  of  wavenumber 
magnitude  and  direction  for  the  generalized  Ekman  shear  in¬ 
stability  problem.  R= 1 1 0 . 0 ,  r*=-0. 0001 . 
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Growth  rate  of  most  unstable  modes  as  a  function  of  normalized 
density  gradient  p*  for  the  generalized  Ekman  shear  instability 
problem.  R=500.0,  ! k ’ =0.5,  e=20°. 
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FIGURE  B-9.  Growth  rate  of  most  unstable  modes  as  a  function  of  normalized 
density  gradient  r*  for  the  generalized  Ekman  shear  instability 
problem  with  Coriofis  terms  included.  R=110.0,  |k=0.5,  c*0  . 
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FIGURE  B-10.  Growth  rate  of  most  unstable  modes  as  a  function  of  scaled 
divergence  for  the  Ekman  shear  instability  problem  with 
coriolis  terms  included.  R= 110.0,  jkj=0.5,  e'O  . 
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WAVENUMBER  MAGNITUDE 


FIGURE  B-ll.  Phase  velocity  as  a  function  of  wavenumber  magnitude  for 

various  modes  corresponding  to  SIGMET-calculated  atmospheric 
parameters  _=0°. 


WAVENUMBER  MAGNITUDE 


FIGURE  B-12.  Growth  rate  of  most  unstable  mode  as  a  function  of  wavenumber 
magnitude  corresponding  to  SIGMET-calculated  atmospheric 
parameters.  e=0  . 


C.  Radar  Properties  of  the  Marine  Atmosphere 

An  apparent  correlation  exists  between  a  surface  radar  duct  and 
certain  convective  phenomena.  Consequently,  it  is  desirable  to  try  to 
bring  convective  instability  of  the  marine  boundary  layer  into  the  same 
theoretical  framework  as  that  of  tropospheric  ducting.  It  is  also  desirable 
to  estimate  the  crossection  for  radar  scattering  from  the  marine  atmosphere. 
Calculations  with  the  SI6MET  1-D  boundary  layer  computer  code  were  described 
in  the  final  report  of  the  first  phase  of  this  work  (Freeman,  1979).  At 
that  time  we  incorporated  a  calculation  of  the  radar  structure  function 
based  on  the  variances  and  covariance  of  temperature  and  water  vapor 
fluctuations  into  SIGMET.  Consequently,  we  were  able  to  associate  profiles 
of  the  radar  structure  function  with  various  environmental  conditions 
governing  the  marine  atmosphere.  An  example  of  this  data  is  shown  in  Fig 
C-l,  taken  from  the  SIGMET  calculation  described  in  Section  A.  The  radar 
structure  function  is  large  at  the  base  of  the  temperature  inversion  and 
in  a  zone  adjacent  to  the  surface. 

The  latter  region  under  appropriate  circumstances  can  result  in 
a  surface  radar  duct  through  channeling  of  the  beam  by  the  gradient  of 
refractivity.  A  quantity  which  measures  the  ability  of  the  atmosphere 
to  support  an  over-the-horizon  radar  duct  is  the  height-modified  microwave 
index  of  refraction 

M  =  n-1  +  |, 

where  Z  is  the  altitude,  and  R  is  the  radius  of  the  earth.  This  quantity 
depends  on  the  index  of  refraction,  n,  which  is  a  function  of  pressure 


p,  temperature  T,  and  humidity  C: 
n  -1  =  y  P  ( 1  +  -y-)  . 

A  negative  gradient  of  the  quantity  M  indicates  that  ducting  will  occur. 
Consequently,  the  altitude  range  through  which  the  gradient  is  negative 
determines  the  depth  of  the  surface  duct.  In  Fig.  C-2  we  show  the  profile 
of  M  early  in  the  SIGMET  calculation  mentioned  above.  This  indicates  that  a 
surface  duct  is  present  having  a  height  of  ~15m,  in  addition  to  an  elevated 
duct  at  abound  500m.  The  M-profile  later  in  the  calculation  (Fig.  C-3)  at 
local  time  =  0108.,  shows  that  the  surface  duct  has  essentially  dissipated 
while  the  elevated  duct  has  risen  in  altitude. 

Recently,  several  calculations  of  the  evolution  of  profiles  of 
microwave  refractive  index  have  been  presented  by  Burk  (1980a),  who  also 
employs  a  1-D  boundary  layer  model.  Several  different  cases  are  examined, 
two  of  which  correspond  to  marine  atmospheres.  These  cases  display  surface 
duct  formation,  although  greater  attention  is  devoted  to  elevated  ducts. 

The  characteristics  of  the  surface  ducts  are  not  resolved  in  the  figures 
presented  by  Burk  and  he  has  not  studied  the  systematics  of  these  ducts. 

We  have  not  had  an  opportunity  to  survey  the  dependence  of  duct 
height  on  atmosphere  parameters  either.  However,  it  is  clear  that  such 
a  study  is  feasible  using  a  tool  like  SIGMET.  In  addition  to  extending 
our  understanding  of  conditions  under  which  ducting  takes  place,  it  could 
be  determined  to  what  degree  the  surface  duct  is  associated  with  the 
susceptabil ity  of  the  atmosphere  to  trigging  of  convection. 

A  somewhat  similar  study  can  be  made  of  the  radar  structure  function. 


5cain,  Burr.  (1530b)  has  recently  given  examples  of  the  profiles  of  the 
crowave  structure  function,  in  addition  to  those  for  the  optical  range 
and  for  acoustics.  Using  the  same  three  cases  for  illustration  as  in 
Burk  (1980a),  the  altitude  dependence  and  relative  contributions  of  temp¬ 
erature  and  moisture  fluctuations  were  evaluated.  Calculations  have  also 
been  made  recently  (Burk,  1980c)  examining  the  effects  of  a  gradient  of 
sea  surface  temperature  on  radar  ducting  and  the  structure  coefficient. 
These  calculations  indicate  that  there  is  a  substantial  change  in  these 
quantities  due  to  air  advection  over  moderate  gradients.  Such  effects 
are  expected  to  be  significant  in  coastal  waters.  Similar  calculations 
were  performed  by  Lewellen  and  Teske  (1975). 

We  can  conclude  (as  in  the  case  of  ducting)  that  theoretical 
models  of  boundary  layer  structure  and  associated  radar  structure  function 
are  available.  The  models  have  not  been  carefully  compared  with  data  or 
evaluated  parametrically.  Investigations  of  the  systematics  of  ducting 
and  scattering  are  desirable  as  part  of  a  program  to  improve  understanding 
of  the  radar  properties  of  the  marine  atmosphere.  In  particular,  radar 
properties  can  be  correlated  with  weather  conditions  affecting  other  local 
phenomena  such  as  the  formation  of  convective  cells. 


nocHMHr® 


> 


I  D.  Internal  Waves  in  the  Marine  Atmosphere 

Under  suitable  environmental  conditions  the  marine  atmosphere  will 
support  the  propagation  of  internal  waves  which  may  be  sensed  in  a  number 

t  of  different  ways.  These  waves  provide  an  important  mechanism  for  the 

transport  of  energy  and  momentum  throughout  the  atmosphere,  they  disturb 
the  mean  flow  in  regions  where  they  are  generated  and  where  they  are  dis- 

*  sipated,  and  they  represent  a  systematic  variability  of  the  atmospheric 
state  which  can  influence  boundary  layer  measurements  (SethuRaman,  et. 
al.,  1980). 

*  It  is  of  interest  to  study  the  properties  of  internal  wave  propaga¬ 
tion  as  a  topic  of  general  interest  in  the  physics  of  the  marine  atomsphere 
in  order  to  assertain  their  rode  in  the  processes  mentioned  above.  More 
specifically,  however,  internal  waves  give  rise  to  atmospheric  disturbances 
which  may  be  detected  by  radar;  we  are  interested  in  examining  whether 
they  could  offer  an  explanation  for  clear  air  echoes  which  propagate  with 

5 

respect  to  the  ambient  wind  field.  It  is  considered  to  be  an  established 
fact  that  large-amplitude  internal  waves  on  breaking  produce  patches  of 
clear  air  turbulence,  and  it  is  likely  that  other  detectable  events  require 

*  large  amplitude  waves  as  well. 

We  are  not  able  to  examine  nonlinear  aspects  of  internal  wave 
propagation  in  this  investigation,  but  the  boundary  layer  linear  perturbation 

1  calculation  described  in  previous  sections  is  useful  for  part  of  such  a 

study.  We  can  take  into  account  the  profile  of  stability,  the  effects  of 
wind  shear,  and  the  damping  by  boundary  layer  turbulence  and  of  convectively 

unstable  regions.  Propagation  speed,  camping  "ate,  and  vertical  distribution 
ox  the  perturbation  can  be  calculated  as  functions  of  wavelength,  propagation 
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direction,  and  mode  number.  It  appears  that  the  internal  wa.e  propagation 
problem  has  not  been  examined  in  this  generality  previously. 


► 

I 

We  have  performed  calculations  -sing  the  PERT  code,  demonstrating 
its  applicability  to  this  class  of  problem.  One  particular  case  of  mea¬ 
surements  of  May  16th,  1979  was  chosen  *or  study;  the  lower  atmosphere  is 
near  neutral  stability  (but  stable)  and  is  capped  by  a  more  stable  region. 

^  The  wind  direction  is  substantially  independent  of  altitude,  and  the  wind 

speed  increases  with  altitude,  rapidly  increasing  near  the  surface  and 
changing  more  gradually  in  the  vicinity  of  the  inversion  base.  The  profiles 
^  of  buoyancy  frequency  and  wind  speed  are  shown  in  Fig.  D-l.  In  this  in¬ 

stance,  the  wind  direction  was  so  nearly  independent  of  altitude  that  the 
turning  with  height  has  been  neglected.  The  buoyancy  freauency  shows 
t  maxima  at  altitudes  well  above  the  base  of  the  inversion;  the  lowest  peak 

is  at  -  900m.  Qualitatively,  for  large  wavenumber  perturbations  the  dis- 
trubance  tends  to  be  concentrated  in  the  vicinity  of  the  buoyancy  peak 
y  and  the  phase  speed  is  small.  Calculations  were  carried  out  for  this 

problem  assuming  that  the  turbulent  dissipation  is  zero.  Wave  directions 
transverse  to  and  aligned  with  the  wind  were  examined.  The  former  cases 
I  are  unaffected  by  the  wind,  and  form  positive/negative  pairs  of  real 

eigenvalues  corresponding  to  waves  traveling  in  opposite  directions. 

When  the  waves  are  aligned  with  the  wind  the  pairs  of  solutions  are  split 
i  apart,  corresponding  roughly  to  the  augmentation  of  the  phase  speed  for  waves 

traveling  down  wind  and  a  reduction  of  phase  speed  for  waves  traveling 
against  the  wind. 

In  Fig  D-2  the  phase  speed  as  e  function  of  the  magnitude  of  the 
wavenumber  of  the  perturbation  is  shown  for  the  transverse  case.  It  is 

» 
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important  to  note  that  the  maximum  phase  speed  is  -  4r/s,  wrich  is  acnieved 
for  small  wavenumber  and  the  first  mode.  Tne  smallest  mode  number  (cor¬ 
responding  to  the  highest  frequency)  will  give  the  largest  degree  of  ducting 
of  the  internal  wave,  resulting  in  localization  of  the  amplitude  around  the 
peak  of  the  buoyancy  frequency.  This  tendency  is  illustrated  in  fig.  D-3 
in  which  the  amplitude  of  the  mode  1  solution  for  k  =  .015  is  shown. 

We  note  that  the  maximum  amplitude  is  found  near  the  height  of  the  peak  of 
buoyancy  frequency,  and  that  the  wave  is  substantially  confined  to  the  region 
of  large  buoyancy  frequency.  For  this  case,  however,  the  phase  speed  is 
no  larger  tham  .7m/s,  which  is  considerably  smaller  than  the  observed 
disturbance  speeds.  With  decreasing  wavenumber  the  phase  speed  increases 
at  the  same  time  that  the  internal  wave  occupies  a  wider  channel.  We  find 
that  the  wave  is  still  substantially  confined  to  the  region  above  the  base 
of  the  inversion  for  a  wavenumber  as  small  as  .00 4  (wavelength  -4500m.), 
for  which  the  phase  speed  is  ~2m/s.  This  speed  is  still  rather  small,  and 
it  is  difficult  to  raise  it  appreciably  for  the  case  in  question.  We  have 
already  chosen  the  fastest  (the  fundanental)  mode  and  we  quote  the  phase 
speed  rather  than  the  group  speed.  Higher  speeds  in  general  are  achieved 
by  a  larger  buoyancy  frequency,  by  longer  wavelengths  of  disturbance, 
and  by  a  deeper  stable  zone  of  the  atmosphere.  For  longer  wavelengths 
for  our  case,  however,  we  tend  to  spread  the  internal  wave  amplitude  throuah- 
out  the  boundary  layer  thickness,  and  it  would  be  necessary  to  invoke  a 
mechanism  (such  as  shear  instability)  for  creating  turbulence  near  the 
inversion  base. 

In  Fig.  D-4  the  corresponding  results  for  phase  speed  vs.  wavenumber 
are  shown  for  the  internal  waves  oriented  parallel  or  anti-parallel  with 
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cne  wind.  Waves  having  1 ow  node  number  propagating  witn  the  wind  are 
o-inarily  rioppler  shifted  to  a  phase  speed  higher  by  acorcximately  4m/s. 

At  higher  wavenumbers  the  effects  of  the  terms  contain-'ng  wind  curvature 
and  slope  become  significant,  resulting  in  a  smaller  decrease  in  phase 
soeed.  For  waves  propagating  against  the  wind  direction  the  dispersion 
curve  shows  larger  departures  from  the  zero-wind  case.  Fig  D-4  indicates 
that,  in  addition  to  a  general  displacement  upwards  by  -  4m/s  of  the  negative 
phase  speed  modes,  there  is  considerable  re-ordering  and  changing  of  shape 
of  the  dispersion  curves.  Some  mode  crossing  cases  occur,  and,  in  contrast 
to  the  transverse  wave  cases,  the  imaginary  part  of  the  freauency  is  not 
zero.  Due  to  the  two-peaked  structure  of  D-l,  it  is  qualitatively  reasonable 
to  expect  some  mode  crossing  due  to  resonance  between  the  two  propagation 
channels.  And  it  is  also  possible  that  there  may  be  shear-unstable  modes 
in  this  case.  However,  the  anti-parallel  modes  almost  invariably  satisfy 
conditions  for  critical  level  occurrence.  It  is  possible,  due  to  the 
coarse  zoning  in  the  vertical,  that  these  cases  void  catastrophic  growth 
(and  absorption),  but  are  nonetheless  inaccurate  (most  of  the  critical  levels 
will  occur  near  the  surface  where  the  amplitude  frequently  is  otherwise 
small).  We  have  not  explored  this  question  in  sufficient  detail  to  evaluate 
the  effect  of  the  critical  levels  on  internal  waves  propagating  against 
the  wind. 

In  order  to  evaluate  the  effect  of  turbulent  dissipation  and  to  take 
a  step  toward  treating  critical  levels  mere  accurately  a  qualitative 
study  of  the  effects  of  the  dissipation  terms  was  made.  We  may  anticipate 
the  trends  resulting  from  dissipation.  For  a  strongly  stable  atmosphere 
dissipation  will  be  localized  near  the  surface  where  wind  shear  produces 
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a  source  of  turbulence.  Internal  waves  having  appreciable  amplitude  at 
low  altitudes  will  be  damped  while  waves  which  are  confined  to  elevated 
regions  of  large  stability  will  be  much  less  affected.  Consequently, 
waves  originating  in  a  distant  source  will  be  filtered  so  that  waves  of 
higher  frequency  predominate.  We  also  expect  that  critical  layer  absorp¬ 
tion  will  affect  the  propagation  of  the  high  frequency  waves  traveling 
in  the  upwind  direction  when  the  profile  of  the  wind  speed  results  in  a 
zone  of  vanishing  intrinsic  frequency.  In  this  case  the  waves  exhibit  fine 
structure  near  the  critical  level  and  enhanced  dissipation  occurs.  We  have 
made  one  calculation  which  accounts  qualitatively  for  dissipation  by  in- 
troducing  a  height-independent  diffusivity  having  a  magnitude  of  2m  /s, 
corresponding  to  a  stably  stratified  atmosphere.  The  dispersion  calculations 
for  aligned  internal  wave  propagation  were  repeated  with  this  addition. 

The  real  phase  speed  vs.  wavenumber  is  shown  in  Fig  D-5;  comparison  with  D-4 
indicates  that  dissipation  has  a  significant  effect  on  phase  speed.  While 
the  lower  modes  propagation  with  the  wind  are  scarcely  affected,  the 
up-wind  modes  are  substantially  changed.  We  also  examined  the  imaginary 
parts  of  the  frequency  and  found  that  no  growing  modes  are  present,  and 
that  significant  wave  decay  takes  place.  For  the  down-wind  modes  the 
decay  is  not  large;  the  first  mode,  for  example,  has  a  3-hour  e-folding 
decay  time  for  k  =  .006  and  longer  lifetime  for  longer  wavelengths.  Higher 
modes  decay  somewhat  faster.  The  up-wind  modes,  on  the  other  hand,  are  strongly 
damped;  the  lowest  phase  speed  mode  for  k  =  .006  decays  in  a  few  minutes, 
although  some  of  the  higher  modes,  paradoxically,  decay  more  slowly,  Although 
we  have  not  explored  this  behavior  systematically,  the  behavior  suggests  that 
the  diffusion  terms  have  the  effect  of  providing  the  damping  required  to  dis- 
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sipate  interna1  wave  energy  at  the  critical  levels. 

Finally,  the  dependence  on  internal  wave  d’ section  .-.as  calcu'eted  for 
the  case  k  =  .001  and  a  diffusivity  of  K  =  2m^/s.  Ir.  ric.  1-6  the  phase 
speed  of  the  first  two  modes  are  plotted  as  a  functior  argle,  ranging 
from  the  transverse  direction  to  the  down  wind  direction.  "he  phase  speeds 
each  shift  by  ~  4m/s  as  expected  from  the  doppler  shi*t.  The  maximum 
phase  speed  of  these  waves  is  not  much  larger  than  4m/ s.  The  effect  of 
wave  direction  is  to  orient  the  waves  with  respect  to  the  wind  so  that 
higher  or  lower  phase  speeds  are  attained.  Turbulent  dissipation  pre¬ 
ferentially  damps  the  higher  mode  and  up-wind  traveling  waves.  The  low 
mode  waves  with  component  of  wave  number  in  the  down-wind  direction  have 
longer  lifetimes;  the  longer  wavelengths  can  survive  many  hours,  but  those 
of  short  wavelength  have  an  e-folding  decay  time  of  -  1  hour. 

The  study  of  internal  wave  propagation  using  the  linear  perturbation 
equations  provides  information  on  the  characteristics  of  waves  favored  by 
environmental  parameters.  For  the  case  investigated  in  our  example  we 
have  calculated  the  dispersion  relation  and  the  vertical  distribution  of  the 
amplitude  of  internal  waves  having  a  range  of  horizontal  wavenumbers. 

The  technique  could  also  be  applied  to  other  profiles  of  wirdspeed  and 
lapse  rate  to  determine  the  dependence  on  these  parameters.  In  particular, 
the  additional  damping  due  to  an  unstable  region  of  the  buoyancy  profile 
could  be  evaluated. 
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FIGURE  D-2.  Internal  wave  phase  speed  as  a  function  of  wavenumber  nagnitude  for 
May  16,  1979  observed  data.  Wave  direction  transverse  to  wind, 
viscosity  =  0.0. 
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Internal  wave  phase  speed  as  a  function  of  wavenumber  magnitude 
for  May  16,  1979  observed  data.  Wave  direction  along  the  wind, 
viscosity  =  0.0. 
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ANGLE  BETWEEN  WAVE  AND  WIND  (deg) 

Interna]  wave  phase  speed  fcr  modes  1  and  2  as  a  function 
of  angle  between  wave  direction  and  wind  direction  for  Kay 
16,  1979  observed  data.  k  =  0.001,  viscosity  =  2.0  m/sec 


Summary  and  Concluding  Remarks 

The  second  phase  of  our  "Investigation  of  Connective  Instability 
in  the  Marine  Boundary  Layer"  has  been  completed.  The  principal  objective 
of  this  study  has  been  to  develop  a  quantitative  boundary  layer  analysis 
technique  requiring  a  small  development  effort  and  having  low  computing 
demands.  Secondary  objectives  have  been  to  consult  with  Naval  Research 
Laboratory  personnel  on  the  design  and  interpretation  of  marine  atmospheric 
measurements,  and  to  improve  understanding  of  marine  convective  phenomena  afe 
sensed  by  radar. 

A  special  report  was  submitted  outlining  considerations  of  a 
theoretical  nature  in  the  conduct  of  marine  atmospheric  experiments  and 
two  visits  to  the  Naval  Research  Laboratory  were  made  to  discuss  current 
problems  with  laboratory  personnel.  The  remaining  topics  are  reported 
in  this  document  where  they  have  been  organized  into  the  four  sections 
discussed  above. 

In  Section  A  the  mean  structure  of  the  boundary  layer  as  displayed 
by  the  SIGMET  computer  code  is  discussed;  in  Section  B,  the  formulation 
and  mathematical  implementation  of  a  computer  code  and  numerical  studies 
of  the  marine  boundary  layer  are  presented.  This  section  contains  the  major 
development  of  the  study.  In  Section  C,  radar  properties  of  the  marine 
atmosphere  are  reviewed  and  in  Section  D  a  study  of  atmospheric  internal 
wave  propagation  is  reported. 

The  tools  resulting  from  this  investigation  and  previous  work  on 
atmospheric  boundary  layer  modeling  form  a  coordinated  set  having  broad 
applications.  The  SIGMET  boundary  layer  computer  code  takes  account  of 
the  processes  affecting  the  evolution  of  the  profiles  of  density,  moisture, 
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Horizontal  wind,  and  diffusivity  in  the  boundary  layer  ;n  response  to 
synoptic-scale  weather  forcing.  In  addition,  quantities  related  to  the 
radar  sensing  of  the  atmosphere  can  be  derived  from  these  profiles. 

Finally,  with  the  development  of  the  linear  perturbation  computer  code 
(PERT),  the  profiles  can  be  employed  as  input  data  to  assess  wave  propagation 
and  stability  of  a  given  boundary  layer.  The  perturbation  code  is  suf¬ 
ficiently  general  to  unify  such  previously  diverse  effects  as  internal 
wave  propagation,  shear  layer  instability,  convective  instability,  critical 
layer  absorption,  and  other  various  boundary  layer  effects  leading  to 
damping  or  instability. 

In  our  investigation  we  have  sampled  a  few  of  the  applications  of 
these  codes.  With  the  boundary  layer  code  the  modified  radar  index  of 
refraction  and  the  radar  structure  function  were  formed.  A  large  number 
of  perturbation  calculations  (~  200)  have  also  been  carried  out  to  test 
several  versions  of  the  PERT  code  and  to  investigate  some  new  aspects  of 
boundary  layer  stability.  We  evaluated  the  effect  of  density  stratification 
on  the  Ekman  unstable  boundary  layer  taking  into  account  the  role  of  the 
Coriolis  term,  and  we  investigated  the  generalized  dispersion  relation 
of  internal  waves  in  a  boundary  layer  containing  a  shearing  wind  field  and 
turbulent  dissipation.  Finally,  we  demonstrated  how  the  mean  boundary- 

layer  code  (SIGMET)  and  the  perturbation  code  can  be  coupled  together 
to  perform  a  boundary  layer  assessment  of  greater  generality. 

Further  use  and  development  of  the  above  tools  was  not  possible 
within  the  scope  of  this  investigation,  but  several  interesting  problems 
can  be  studied  in  the  future.  Particularly  valuable  is  the  possibility 
of  correlating  several  different  effects  through  the  combined  use  of  these 
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codes.  We  can  determine,  for  example,  the  atmospheric  conditions  under 
which  radar  surface  ducting  takes  place  and  the  associated  boundary  layer 
stability.  Another  problem  of  interest  is  to  determine  the  internal  wave 
properties  of  the  boundary  layer  for  a  range  of  weather  conditions. 

Clearly,  a  linear  stability  analysis  reveals  only  one  aspect  of  the 
secondary  circulations  of  the  boundary  layer  which  result  from  instability. 
To  understand  the  state  of  the  atmosphere  reached  after  these  instabilities 
have  grown  requires  that  nonlinear  effects  be  included.  Since  conventional 
methods  of  nonlinear  fluid  dynamics  are  cumbersome  and  very  expensive, 
it  is  desirable  to  retain  the  simplicity  and  generality  of  the  linear 
analysis  in  a  nonlinear  investigation.  Some  hope  for  accomplishing  this 
goal  is  held  out  by  expansion  methods  restricted  to  a  few  terms. 
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